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' We discuss phase coexistence of polydisperse colloidal suspensions in the presence of adhesion forces. The combined effect of polydispersity 
and Baxter's sticky-hard-sphere (SHS) potential, describing hard spheres interacting via strong and very short-ranged attractive forces, 
give rise, within the Percus-Yevick (PY) approximation, to a system of coupled quadratic equations which, in general, cannot be solved 
either analytically or numerically. We review and compare two recent alternative proposals, which we have attempted to by-pass this 
, difficulty. In the first one, truncating the density expansion of the direct correlation functions, we have considered approximations 
, simpler than the PY one. These Cn approximations can be systematically improved. We have been able to provide a complete analytical 
description of polydisperse SHS fluids by using the simplest two orders Cq and Ci , respectively. Such a simplification comes at the price of 
■ a lower accuracy in the phase diagram, but has the advantage of providing an analytical description of various new phenomena associated 
' with the onset of polydispersity in phase equilibria (e.g. fractionation). The second approach is based on a perturbative expansion of 
' the polydisperse PY solution around its monodisperse counterpart. This approach provides a sound approximation to the real phase 
behavior, at the cost of considering only weak polydispersity. Although a final seattlement on the soundness of the latter method would 
require numerical simulations for the polydisperse Baxter model, we argue that this approach is expected to keep correctly into account 
, the effects of polydispersity, at least qualitatively. 

1 Introduction 

' New technological advances in physico-chemical manipulation of colloidal mixtures have brought up again 
the issue of theoretically understanding the phase behaviour of polydisperse systems (jl|). 'Polydispersity' 
in colloidal solutions means that, due to their production process, suspended macroparticles with the same 
chemical composition cannot be exactly identical to each other, but in general have different sizes, and 
possibly different surface charges, shapes, etc. In practice, a polydisperse system can be reckoned as a 
mixture with very large - or essentially infinite - number M of different species or components, identified 
by one or several parameters (M large but finite refers to discrete polydispersity, whereas M ^ oo with a 

! continuous distribution of polydisperse parameters corresponds to continuous polydispersity). The present 

' paper will consider discrete polydispersity of spherical colloidal particles, with their diameter being the 

. only polydisperse attribute (size-polydispersity). 

' When polydispersity is not negligible, the phase behaviour becomes much richer, but the determination 
of phase transition boundaries requires a much more involved formalism, compared to the monodisperse 
counterpart. In fact, the coexistence condition in terms of intensive variables requires that all phases 
must have equal temperature, pressure and chemical potentials of the M components. In the presence of 
polydispersity, one should thus solve a number of equations of the order of M^, a task which is practically 
impossible for M large or infinite. 

However, the study of phase equilibria can conveniently start from the appropriate thermodynamic po- 
tential, which is the Helmholtz free energy A when the experimentally controlled variables are temperature, 
volume and numbers of different colloidal species. In the one-component case, the coexistence condition of 
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equal pressure and chemical potential has a simple geometrical interpretation in terms of free energy den- 
sity a: the densities of two coexisting phases are determined by constructing a double-tangent to a plotted 
versus particle density. This recipe leads to the well known Maxwell construction, which connects suitably 
selected points along a van der Waals (vdW) subcritical isotherm, in order to 'reduce' its unphysical loop 
to a constant-pressure line characteristic of a first-order phase transition. 

In the polydisperse case, a significant progress in the very difficult problem of predicting phase equilibria 
can be obtained for models with truncatable free energies (jll). Here 'truncatable' means that the excess 
free energy of the polydisperse system turns out to depend only on a finite number of moments of the 
distribution corresponding to the polydisperse attribute (the diameter a in the simplest cases). For spherical 
colloids, the excess free energy of the vdW model extended to polydisperse fluids has such a truncatable 
structure. Due to this property, this vdW theory has often been employed as the simplest model to 
investigate the effects of polydispersity on the gas-liquid transition On the other hand, the influence 

of polydispersity on freezing has been addressed by using the hard-sphere (HS) mixture model, which 
also admits a truncatable free energy (1) (for the fluid phase, the Boublik-Mansoori-Carnahan-Starling- 
Leland (BMCSL) (js') expression was employed). It is worth recalling that it is currently believed that size- 
poly dispersity might destabilize crystallization, eventually inhibiting freezing above a certain 'terminal' 
value of polydispersity ([l|). 

The present paper focuses on - and reviews ~ a number of recent attempts to investigate polydisperse 
phase equilibria, at least within some approximations, for another prototype model useful for studying 
colloidal suspensions, namely Baxter's sticky-hard-spheres (SHS) (4). Here the particles are hard spheres 
with a surface adhesion, and the corresponding potential can be obtained as a limit of an attractive square- 
well which becomes infinitely deep and narrow, according to a particular prescription which ensures a finite 
non-zero contribution of adhesion to the second virial coefficient ('sticky limit'). For the one-component 
version of this model, Baxter and collaborators Q) solved the Ornstein-Zernike (OZ) integral equation 
coupled with the Percus-Yevick (PY) approximation ('closure'). This fully analytical solution allows to 
determine all structural and thermodynamical properties of the SHS fluid. On the other hand, the multi- 
component PY solution, which soon followed Baxter's work ([1;@), is practically inapplicable in the presence 
of significant polydispersity. In fact, it requires the computation of a set of parameters {A^j} determined 
by a system of M(M + l)/2 quadratic equations, which - in general - cannot be solved even numerically 
for a mixture with a large number of components. Moreover, even in the presence of a general solution for 
this non-linear algebraic system, the problem of phase coexistence would still remain out of reach in view 
of the previous remarks. 

In a series of recent papers (0; H; S [HI; fisl ). we have attempted to make some progress along two 
different lines. 

First, starting from the density expansion of the cavity function at contact, we have considered a sequence 
of simpler approximations (compared to the PY one) fn: Isl: lil: [lol: [ill) . Within the two simplest ones among 
these approximations, denoted as Cq and Ci (for reasons which will become clear in the following), we have 
been able to derive analytically all relevant information regarding structure and thermodynamics, including 
the phase coexistence, in view of the fact that the corresponding free energy turns out to be truncatable 



(jlll ). Due to the simplicity of Cq and Ci, it is however reasonable to expect these approximations to fail 
at high packing fractions, with a consequently incomplete or even incorrect description of the effects of 
polydispersity on the phase diagram. 

Therefore, in collaboration with Peter Sollich, we have recently explored a second approach (|l3), where 
the expansion variable (which must be small) is an appropriate polydispersity index. In such a way, we have 
tried to solve the non-linear algebraic system - involved in the PY result - perturbatively in polydispersity, 
starting from the monodisperse PY solution. 



2 Baxter's SHS model and PY solution 



The SHS model is defined as limiting case of a particular square- well (SW) model (3), based upon a 
potential including steeply repulsive core and short-ranged attractive tail, i.e. 
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where cTj is the HS diameter of species i, ef- 
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where the conventional Baxter parameter Tij is an unspecified increasing function of T, while r^~^ measures 
the strength of surface adhesion or 'stickiness' between particles of species i and j. 

The 'sticky limit' of ^^.^xtsf SW ^^-j corresponds to taking Wij — > 0. While the SW width goes to zero, its 
depth e^'^^*'^'' diverges, giving rise to a Dirac delta function in the Boltzmann factor 



1. e. 



-mri-) = 9{r- a,j) + tij a^,5 (r - a,j) 



(4) 



where P = {kBT)~^ , while 6 and 5 are the Heaviside step function and the Dirac delta function, respectively. 

The advantage of the sticky limit is that one effectively deals with a single parameter Tij for each 
pair, rather than a combination of energy and length scales (as occurs in the square-well model, for 
which no analytical solution is known). On the one hand, this particular limit has the disadvantage of 
introducing some pathologies in the model, notably in the one-component case On the other hand, 
Baxter's model represents the simplest paradigmatic example of a combination of steep repulsion and 
short-range attraction which entails a complete analytical solution in the one-component case, within a 
robust approximation such as the PY closure. 

In the multicomponent case, the PY solution of the OZ equation in terms of Baxter's factor correlation 
function reads (j^; M) 



I 0, elsewhere 
where the expressions for the parameters Oj and bi may be found in (0), while the quantity 



^ij — tij yij{(yij) O'ij , 



(6) 



which depends on the cavity function at contact yij{(Jij), must be solution of the following system of 
quadratic equations 



A, 



a 



+ A. E 



1,2, 



(7) 



Here, Xm is the molar fraction of the m-th species (m = 1,...,M), while Oij = tij y^^~^^ (cij) crfj, 
0ii = 12/9 tijaij (p is the total number density), and Tij = aiaj/{l — ij), with r] being the packing fraction 
(|l2l ). The solution of these equations for {Aj^} is the real bottleneck of the multi-component PY result, 
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as mentioned in the Introduction: for large M (and in particular for M ^ oo ) this calculation is next to 
impossible, neither analytically nor numerically. 

As a consequence, although the PY closure is commonly believed to be very sound for short-range 
potentials (for one-component SHS fluids this was confirmed by recent numerical simulations (15)), one 
has to conclude that, in the multi-component (polydisperse) case, the PY solution has a very limited 
practical usefulness, since its solution scheme cannot be fully accomplished. This is the reason why other 
possible routes have been attempted, as we discuss next. 



3 Simplified closures: the class of Cn approximations 

A 'closure' is a relationship, added to the OZ equation, between the direct correlation function Cij{r) and 
hij{r) = gij{r) — 1 or ^ij{r) = hij{r) — Cij{r) {gij{r) being the radial distribution function) (jiol ). 

Let us go back to Baxter's SW model given by Eq. ([1]) (i.e. before the 'sticky limit'), and consider the 
following general class of 'closures' (0) 



-1 - (r) < r < cjjj 

^shi' 



Ci,{r) = { c"(r) a,,<r<R,^ (8) 



r > Rij . 



The expression Cij{r) — —1 — (r) inside the core (r < Uij) is exact and dictated by the HS potential. 
The form outside the well (r > Rij) may then be identified with the PY approximation, 



cl^{r)=Mr)[l + l^,{r)], (9) 

since, for Baxter's potential, the Mayer function, fij{r) = exp [—(3(j)ij{r)] — 1, vanishes for r > Rij. 

The choice of c-^'^™'^(r) inside the well (region which 'shrinks' in the sticky limit) defines one particular 
closure within the proposed class. Of course, cfj™'^(r) = cf^{r) corresponds to the PY approximation. On 
the other hand, when c^hrmk^- _^ cfj^{r), we are in the presence of mixed closures, which have frequently 
appeared in the literature (|l9l ). In order to define mixed closures simpler than the PY approximation, we 
consider the density expansion of the exact direct correlation function (16), and denote as C„ approxima- 
tion a truncation of this series to order 0(/9"). The simplest two approximations are 



dtf''Hr) = Mr) (Co closure) 

^shrink(^) = (r) [1 + ( Efc Pk hk * fko) {r)] (Ci closure), ^ ^ 

where pk is the number density of species k, while * denotes convolutive integration |iqI). 

In the 'sticky limit' Rij — >• af- the well region shrinks, but a 'memory' of the approximation chosen for 
^shrink j-gj-Q^ins in the solution of the OZ integral equation. In fact, although all solutions qij{r) correspond- 
ing to closures belonging to the class given by Eq. ([8]) have the same functional form as the PY solution 
- Eq. dS]) " , each closure is characterized by its own approximation to yij{crij), which is involved in the 
expressions of the parameters Oj, bi, Aij. For instance, the Co and Ci approximations correspond to 



Vijic^ij) = 1 (Co closure) 

Viji^ij) = 1 + yifi(^ij) V (Ci closure), 

which are, respectively, the zeroth- and first-order truncations of the density expansion for the exact 
cavity function at contact (see Ref. (jiol ) for details). 

While a brute-force truncation of the above-mentioned density expansions leads to analytical expressions 
simple enough to be applied to the multi-component (polydisperse) case, one should reasonably expect 
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less accuracy, expecially in the high-density regime. In the one-component case, we can carefully check 
this point. 

In Figure [T] coexistence curves obtained from the Cq and Ci approximations are compared with the PY 
ones (using both compressibility and energy routes), and with Monte Carlo simulations from Ref. (IS). It 
is apparent how the PY energy route (PYE) yields a rather precise representation of the MC results, unlike 
the compressibility route (PYC). It is worth noting that the results stemming from the Ci approximation, 
although rather close to the MC data in the low-density branch, clearly fail to accurately reproduce them 
for higher densities, as expected. 

In spite of their lack of accuracy, the Cq and Ci approximations provide however a rather sound basis 
for getting some insight into phase equilibria of polydisperse SHS fluids, since they allow simple analytical, 
or semi-analytical, treatments. 

A first important feature of the Cq and Ci approximations for polydisperse SHS is that the corresponding 
free energy has a truncatable structure, that is it depends upon few (4 at most) moments of the (discrete) 
size distribution, = {'^/6)pJ2j^j^j with 1^ = 0, 1,2,3. 

A second remarkable fact is that the Cq and Ci approximations are able to describe the so-called 
fractionation phenomena characteristic in phase equilibria of polydisperse systems. While we refer to a 
recent review (|l|) for detailed description of the increased complexity in the polydisperse phase diagrams, 
here we just mention the two important points. First, fractionation means that daughter phases, obtained 
from demixing of a parent homogeneous phase, need not have the same composition of the parent phase. As 
a consequence, there is no a single coexistence line ('binodal') as in the one-component case, but one rather 
finds a cloud curve, representing the temperature-density dependence line of the low-density majority phase 
('gas'), and a shadow curve representing the temperature-density dependence of the high-density minority 
phase (incipient 'liquid'). While for one-component systems these two curves are identical, for polydisperse 
systems in general they are not, with the exception of the critical point where they coincide by definition. 

However, in order to apply the Cq and Ci approximations to the multi-component SHS model, we have 
to tackle a further important problem, that is the definition of the stickiness parameters r^. 



4 Size-dependence of stickiness parameters 

In mixtures, Tij will depend on the particular pair i,j considered, and reasonably should be expected to 
be related to the particles sizes. Assuming that we are dealing only with size-polydispersity, we can always 
decouple temperature and adhesion as 



1 



1 



1 



(12) 



where the last equality stems from the assumption of size-polydispersity and of a purely pairwise potential. 
Unfortunately, the exact form of the size-dependence of these stickiness parameters is still an open problem, 
due to the lack of experimental and theoretical insights on this (fl3l ). On the other hand, few guidelines 



- based on arguments discussed in Refs. (jll 
dependences 



12l ) - provide, as reasonable and plausible, the following 



' o-'^/afj Case I , 
aiaj/afj Case II , 
1 Case IV 

, co/cTij Case V . 



(13) 



Here <to is a characteristic reference length (e.g. the parental mean diameter) and the numbering of the 
various cases follows the convention of previous work (111; 1^ . 

Figure [2] reports the results of the calculation of the cloud and shadow curves for polydisperse SHS 
within the simple Cq approximation. Here and below, polydispersity is measured by an index s, which is 
the normalized standard deviation of the size distribution. Hence s = corresponds to a mono-disperse 
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case, whereas s = 0.1 and s = 0.3 indicate moderate and significant polydispersity, respectively. The top 
panel of Fig. [2] depicts the results for case I of size-dependence of stickiness parameters. As s increases, 
the coexistence region shrinks, thus suggesting that polydispersity disfavours the phase transition. On the 
other hand, this trend is markedly case-dependent, as illustrated in the bottom panel of Fig. [21 where 
the cloud-shadow pairs with polydispersity s = 0.3 are displayed for different cases of size-dependence. It 
can be clearly seen that while, for cases I and V the same trend is observed, case IV seems to suggest a 
widening of the phase coexistence region (and hence a favouring of the phase transition). 

In view of the lack of numerical simulation for polydisperse SHS to compare with, we have no way, at 
the present stage, to check how realistic those results are. On the other hand, we might suspect, based on 
the comparison in the one-component case, Cq to fail to provide accurate representation in the region of 
low temperatures and high densities. This is the reason why other possible approaches have been recently 
tested. We now illustrate a different perturbative approach, which has proved to be promising in this 
context. 



5 Perturbative treatment of the SHS-PY solution 



The main difficulty in dealing with the FY solution for polydisperse SHS stems from the solution of the 
coupled quadratic system of equations ([7]), as remarked. As the one-component case has a well-defined 
solution, one might then suspect that - for weak polydispersity ~ a perturbative expansion around this 
solution might include the main effects of polydispersity. This is in fact what happens, as recently shown 
( 12) by exploiting a general perturbation theory due to Evans (flTl). The main idea is that, for weak 
polydispersity, size-distributions are narrowly peaked around a mean reference value (do in the present 
case), and hence all quantities such as 



= « 1, (14) 

o"o 

are small. Therefore one might expand both eij, and all quantities appearing in Ajj, in powers of 5i. Similar 
expansion can be performed in the free energy, and hence all thermodynamic quantities can be computed. 
The entire procedure is described in detail in Refs. (fll: [l3: . Here we just summarize the main results. 

The approximate range of validity of the perturbation expansion can be envisaged by considering the 
polydisperse HS case where the 'exact' BMCSL approximation (3) can be compared with the corresponding 
perturbative solution based on the one-component (s = 0) counterpart. This is reported in the top-left 
part of Figure [31 where the quantity l3Pvo (vq = 7rcrg/6) is plotted against the packing fraction ry for 
increasing values of polydispersity. It is apparent how the perturbative solution remains rather close to 
the 'exact' polydisperse BMCSL solution even for moderate polydispersity s < 0.3, which is the maximum 
value considered in the remaining part of the work. The remaining plots in Figure [3] display the effect of 
polydispersity on the FY pressure equation of state as obtained from the energy route and for decreasing 
values of the temperature r. In the one component case s = 0, a van der Waals loop starts to appear 
when we cross the critical temperature Tc ~ 0.1185 coming from the high r regime. Obviously this signals 
the onset of a liquid-gas phase transition, and the corresponding phase diagram can be obtained by a 
standard Maxwell construction by connecting appropriate points with the same pressure. In the presence 
of polydispersity (here represented by choice IV for size-dependence of the stickiness parameters) the same 
procedure cannot be applied due to fractionation, as already discussed. Nevertheless, we can clearly see 
that as s increases, the van der Waals loop region (when present) expands, thus suggesting that phase 
transition is favored by the presence of polydispersity. A similar feature occurs for the polydisperse van 
der Waals model (j3) and for the numerical results of the FY compressibility equation of state (|l^ ) (note 
that in the latter a gap rather that a loop is signalling the onset of the transition) . A somewhat surprising 
feature is that, at fixed packing fraction rj, the pressure decreases with increasing polydispersity less in the 
presence of adhesion rather than in its absence (i.e. for the HS case). An intuitive plausible interpretation 
of this feature can be found in Ref . (jH) . 

The same perturbative approach allows the determination of the cloud and shadow curves for the various 
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cases of size-dependence of r^^^. This is reported in Figure H] for cases II, IV (top panel) and I, V (bottom 
panel) . In the first case the cloud and shadow lines collapse into a single curve, and this can be understood 
on the basis of the particular scaling properties of the free energy to this order in perturbation theory 
In all cases there is a breakdown of the perturbation theory on approaching the critical point, and this is 
a known general drawback of Evans' perturbative scheme. Nevertheless, in all cases and to this order in 
perturbation theory, there is a tendency of the phase coexistence region to increase with polydispersity, 
in qualitative agreement with the intuitive picture obtained from Figure [3l 

It is worth stressing the difference with previous non-perturbative results stemming from the Co solution, 
where all different cases (with the notable exception of IV) were predicting a reduction of the phase 
coexistence region. While in the Co description we have provided a careful treatment of polydispersity 
at the expenses of accuracy of the exploited approximation, in the perturbative description of the PY 
solution, polydispersity is assumed to be small and hence one might suspect that solutions with large 
polydispersities cannot fit within this picture. On a balance, nevertheless, we would favour the latter 
rather than the former description. An almost correct representation of the one-component counterpart, 
is a necessary requirement to check the effect of polydispersity on it, and we are not aware of any physical 
or experimental system where the effects of polydispersity are so strong that they could not kept into 
account, at least at the simplest qualitative level, by the perturbative scheme proposed here. Along this 
line, some further proposals have been put forward in Ref. ((l3) to derive a phenomenological BMCSL-like 
approximation for SHS, which might be regarded as our 'best and simplest guess' to the exact phase 
behavior of polydisperse SHS. Even on the size-dependence of t^^ some possible support of the proposed 
forms may be argued (H; 13). 



6 Conclusions 



In this work we have summarized recent advances in predicting theoretically the phase diagram for poly- 
disperse suspensions of colloidal particles with surface adhesion, within the simple description of Baxter's 
model. Emphasis was put on the crucial - unsolved - step required to get the multi-component SHS-PY 
solution, and the proposed recipies to deal with this problem. This first one is based on a simplification of 
the closure. It has the advantage of allowing a complete analytical analysis on the effects of polydispersity, 
including fractionation, but has the disadvantage of a very questionable accuracy. The second one is based 
on a perturbative method, starting from the energy PY one-component solution, which is known to pro- 
vide an accurate description of the phase diagram. The drawback of this scheme is that it works for mild 
polydispersity, and that cannot describe the changes of the critical point region. Notwithstanding these 
limitations, this novel approach is expected to find practical application in the interpretation of all those 
phenomena where Baxter's model and polydispersity both play a privileged role. 
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Figure 1. Coexistence (binodal) curves for the one-component Baxter model. Both compressibihty (PYC) and energy (PYE) equation 

of state as obtained from the Percus-Yevick approximation (see Ref iQ)) are reported and compared with Co and Ci approximations 
from Ref. llllh and with Monte Carlo simulation (MC) from Ref. In the MC case the continuous line is simply a guide to the eye. 
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Figure 2. (Top) Cloud and shadow curve for model I within the Co approximation at increasing polydispersity: s = 0, s = 0.1 and 
s = 0.3 . (Bottom) Same as above for fixed value of polydispersity s = 0.3 and different choice of the stickiness adhesion (model I, IV 

and V) 
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Figure 3. Behavior of the energy equation of state within our pcrturbative scheme. In all cases the quantity /3Pvo is plotted against 
the packing fraction rj. In clockwise order, the first curve (left, top) reports a comparison of the perturbative versus the 'exact' BMCSL 
solution in the equation of state for polydisperse HS (t = +oo). The other curves report the perturbative solution for the energy 

equation of state within the PY approximation for SHS Baxter model. Results are depicted for three values of temperature 
r = 0.15 > Tc, T = 0.1185 ~ Tc and r = 0.1 < Tc and for different degrees of polydispersity. The choice for the size-dependence of 

stickiness parameters corresponds to model IV. 
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Figure 4. Cloud/Shadow pairs from the pcrturbative results for the PY solution of SHS Baxter model. In clockwise order the results 
for choices II, IV (top) and I,V (bottom) are depicted. In the two top panels, the cloud and shadow curves coincide to this order in 
perturbation, whereas in the bottom panels they axe different. In order to have all pictures on the same scale, the selected value for 

polydispersity is s = 0.3 for models II, IV (top) and s = 0.1 for models I, V (bottom). In all cases the continuous curve represents the 

monodisperse (s = 0) result. 



